#=
产生相应的密度矩阵元
=#


"""
某一个时刻的密度矩阵
"""
function spinful_rho_at(upbds1::Matrix{T}, dnbds1::Matrix{T},
    upbds2::Matrix{T}, dnbds2::Matrix{T}, nebds::Matrix{T}, volume::Vector{Float64};
    λe=_ThermWL, μee=_WLElec, μep=_WLNeur, Cne=_NeurCoul, Δβ=_Δβ) where {T}
    #ρ=ρ^U ρ^K
    #势能求和
    #ee部分
    rpos1 = vcat(upbds1, dnbds1)
    ϕee1 = Kelbg_repulse_pairs(rpos1; λ11=μee)
    rpos2 = vcat(upbds2, dnbds2)
    ϕee2 = Kelbg_repulse_pairs(rpos2; λ11=μee)
    potee = 0.5ϕee1 + 0.5ϕee2
    #ep部分
    ϕep1 = Kelbg_attract_pairs(rpos1, nebds; λ22=μep, Cne=Cne)
    ϕep2 = Kelbg_attract_pairs(rpos2, nebds; λ22=μep, Cne=Cne)
    potep = 0.5ϕep1 + 0.5ϕep2
    pot = potee + potep
    #K部分
    kinup = free_kinetic(upbds1, upbds2, volume; λe=λe)
    kindn = free_kinetic(dnbds1, dnbds2, volume; λe=λe)
    kin = -log(kinup) -log(kindn)
    return Δβ*pot, kin
end


"""
最后一个时刻的交换矩阵
"""
function spinful_rho_exch(upNt::Matrix{T}, dnNt::Matrix{T},
    up0::Matrix{T}, dn0::Matrix{T}, nebds::Matrix{T}, volume::Vector{Float64};
    λe=_ThermWL, μee=_WLElec, μep=_WLNeur, Cne=_NeurCoul, Δβ=_Δβ) where {T}
    #
    rpos1 = vcat(upNt, dnNt)
    ϕee1 = Kelbg_repulse_pairs(rpos1; λ11=μee)
    rpos2 = vcat(up0, dn0)
    ϕee2 = Kelbg_repulse_pairs(rpos2; λ11=μee)
    potee = 0.5ϕee1 + 0.5ϕee2
    #ep部分
    ϕep1 = Kelbg_attract_pairs(rpos1, nebds; λ22=μep, Cne=Cne)
    ϕep2 = Kelbg_attract_pairs(rpos2, nebds; λ22=μep, Cne=Cne)
    potep = 0.5ϕep1 + 0.5ϕep2
    pot = potee + potep
    #交换部分
    emup = exch_free_kinetic(upNt, up0, volume; λe=λe)
    emdn = exch_free_kinetic(dnNt, dn0, volume; λe=λe)
    absdetup = abs(det(emup))
    absdetdn = abs(det(emdn))
    return Δβ*pot, -log(absdetup) -log(absdetdn)
end


"""
某条路径产生的权重
"""
function spinful_rho_path(fbup::WrdFb{T}, fbdn::WrdFb{T}, fbne::WrdFb{T}, volume::Vector{Float64};
    λe=_ThermWL, μee=_WLElec, μep=_WLNeur, Cne=_NeurCoul) where {T}
    Np = length(fbup)# dn 
    Nt, Nd = size(fbup[1])# dn
    #
    rhov = zeros(T, Nt)
    rhop = zeros(T, Nt)
    rhok = zeros(T, Nt)
    for it in Base.OneTo(Nt-1)
        upbds1 = beads(fbup, it)
        upbds2 = beads(fbup, it+1)
        dnbds1 = beads(fbdn, it)
        dnbds2 = beads(fbdn, it+1)
        #电子之间的势能
        #ϕ -> 0.0
        #for i<j in up,dn -> ϕ += 1/2 {ϕ_Kelbg(ri,rj) + ϕ_Kelbg(r'i,r'j)}
        #或者 for i,j in up,dn -> ϕ += 1/4 {ϕ_Kelbg(ri,rj) + ϕ_Kelbg(r'i,r'j)}
        #因为每个都要计算两次，也可以是
        #for i,j in up,dn -> ϕ += 1/2 ϕ_Kelbg(ri,rj)
        rpos1 = vcat(upbds1, dnbds1)
        ϕee1 = Kelbg_repulse_pairs(rpos1; λ11=μee)
        rhov[it] = ϕee1
        #电子质子之间的势能
        #ϕ -> 0.0
        #for i in up,dn for j in ne -> ϕ += 1/2 {ϕ_Kelbg(ri,rj) + ϕ_Kelbg(r'i,r'j)}
        nebds1 = beads(fbne, it)
        nebds2 = beads(fbne, it+1)
        ϕep1 = Kelbg_attract_pairs(rpos1, nebds1; λ22=μep, Cne=Cne)
        rhop[it] = ϕep1
        #电子的动能
        kinup = free_kinetic(upbds1, upbds2, volume; λe=λe)
        kindn = free_kinetic(dnbds1, dnbds2, volume; λe=λe)
        rhok[it] = log(kinup) + log(kindn)
    end
    #最后一个会转回来
    upbds1 = beads(fbup, Nt)
    upbds2 = beads(fbup, 1)
    dnbds1 = beads(fbdn, Nt)
    dnbds2 = beads(fbdn, 1)
    rpos1 = vcat(upbds1, dnbds1)
    ϕee1 = Kelbg_repulse_pairs(rpos1; λ11=μee)
    rhov[Nt] = ϕee1
    nebds1 = beads(fbne, Nt)
    nebds2 = beads(fbne, 1)
    ϕep1 = Kelbg_attract_pairs(rpos1, nebds1; λ22=μep, Cne=Cne)
    rhop[Nt] = ϕep1
    #这个时候就需要交换矩阵了
    #kinup = free_kinetic(upbds1, upbds2; λe=λe, vol=volume)
    #kindn = free_kinetic(dnbds1, dnbds2; λe=λe, vol=volume)
    exchup = exch_free_kinetic(upbds1, upbds2, volume; λe=λe)
    exchdn = exch_free_kinetic(dnbds1, dnbds2, volume; λe=λe)
    absexup = abs(det(exchup))
    absexdn = abs(det(exchdn))
    rhok[Nt] = log(absexup) + log(absexdn)
    #
    return rhov, rhop, rhok
end


"""
将某个位置的粒子移动
"""
function move_at(aip::Int, ait::Int, fb::WrdFb{T}, volume::Vector{Float64};
    λe=_ThermWL) where {T}
    #每次修改的时候，影响it和it-1两个位置的数值
    #<R1|e^h|R2><R2|e^h|><|...<RNt|e^h|R1>
    Np = length(fb)
    Nt, Nd = size(fb[1])
    wfb = WrdFb(T, Np, Nt)
    for ip in Base.OneTo(Np); for it in 1:1:Nt
        wfb[ip][it, :] = fb[ip][it, :]
    end; end
    #
    off = (rand(3).-0.5)*λe
    #println(aip, ait)
    wfb[aip][ait, :] = map_to_1st(wfb[aip][ait, :]+off, volume)
    return wfb
end


"""
对电子的路径进行更新
"""
function update!(fbup::WrdFb{T}, fbdn::WrdFb{T}, fbne::WrdFb{T}, volume::Vector{Float64};
    λe=_ThermWL, μee=_WLElec, μep=_WLNeur, Cne=_NeurCoul) where {T}
    Ne = length(fbup)#dn dn
    Nt, Nd = size(fbup[1])#dn
    Nc = length(fbne)
    #更新位置，第一个和最后一个的更新会涉及交换矩阵，需特殊处理
    for it in Base.OneTo(Nt)
        nebds1 = beads(fbne, it)
        #
        lupt = beads(fbup, it == 1 ? Nt : it-1)
        nupt = beads(fbup, it)
        rupt = beads(fbup, it == Nt ? 1 : it+1)
        ldnt = beads(fbdn, it == 1 ? Nt : it-1)
        ndnt = beads(fbdn, it)
        rdnt = beads(fbdn, it == Nt ? 1 : it+1)
        #println("beads before $(nupt)")# dn
        #dn 
        #println("update it: $(it) lit: $(it==1 ? Nt : it-1) rit: $(it==Nt ? 1 : it+1)")
        #更新前的密度矩阵
        if it != 1
            #<t-1| ρ |t>
            lpot, lkin = spinful_rho_at(lupt, ldnt, nupt, ndnt, nebds1, volume;
            λe=λe, μee=μee, μep=μep, Cne=Cne)
            #println("ρ($(it-1)): pot=$(lpot), kin=$(lkin)")
        else
            #<Nt| ρ |Perm(1)> 
            lpot, lkin = spinful_rho_exch(lupt, ldnt, nupt, ndnt, nebds1, volume;
            λe=λe, μee=μee, μep=μep, Cne=Cne)
            #println("ρ̃($(it-1)): pot=$(lpot), kin=$(lkin)")
            #lpota, lkina = spinful_rho_at(lupt, ldnt, nupt, ndnt, nebds1, volume;
            #λe=λe, μee=μee, μep=μep, Cne=Cne)
            #println("ρ($(it-1)): pota=$(lpota), kina=$(lkina)")
        end
        if it != Nt
            # <t| ρ |t+1>
            rpot, rkin = spinful_rho_at(nupt, ndnt, rupt, rdnt, nebds1, volume;
            λe=λe, μee=μee, μep=μep, Cne=Cne)
            #println("ρ($(it)): pot=$(rpot), kin=$(rkin)")
        else
            # <Nt| ρ |Perm(1)>
            rpot, rkin = spinful_rho_exch(nupt, ndnt, rupt, rdnt, nebds1, volume;
            λe=λe, μee=μee, μep=μep, Cne=Cne)
            #println("ρ̃($(it)): pot=$(rpot), kin=$(rkin)")
            #rpota, rkina = spinful_rho_at(nupt, ndnt, rupt, rdnt, nebds1, volume;
            #λe=λe, μee=μee, μep=μep, Cne=Cne)
            #println("ρ($(it)): pot=$(rpota), kin=$(rkina)")
        end
        #尝试更新位置 up
        for ip in Base.OneTo(Ne)
            _nupt = copy(nupt)
            ran = rand(3).-0.5
            ran = ran*0.1
            #println(ran)
            _nupt[ip, :] = map_to_1st(_nupt[ip, :] + (ran)*λe, volume)
            #println(_nupt, " λe = $(λe)")#dn
            #println("update site $(ip) up: $(nupt) -> $(_nupt)")#dn dn dn dn
            #更新后的密度矩阵
            if it != 1
                #<t-1| ρ |t>
                _lpot, _lkin = spinful_rho_at(lupt, ldnt, _nupt, ndnt, nebds1, volume;
                λe=λe, μee=μee, μep=μep, Cne=Cne)
                #println("ρ($(it-1)): pot=$(_lpot), kin=$(_lkin)")
            else
                #<Nt| ρ |Perm(1)> 
                _lpot, _lkin = spinful_rho_exch(lupt, ldnt, _nupt, ndnt, nebds1, volume;
                λe=λe, μee=μee, μep=μep, Cne=Cne)
                #println("ρ̃($(it-1)): pot=$(_lpot), kin=$(_lkin)")
                #lpota, lkina = spinful_rho_at(lupt, ldnt, _nupt, ndnt, nebds1, volume;
                #λe=λe, μee=μee, μep=μep, Cne=Cne)
                #println("ρ($(it-1)): pota=$(lpota), kina=$(lkina)")
            end
            #右侧的密度矩阵
            if it != Nt
                # <t| ρ |t+1>
                _rpot, _rkin = spinful_rho_at(_nupt, ndnt, rupt, rdnt, nebds1, volume;
                λe=λe, μee=μee, μep=μep, Cne=Cne)
                #println("ρ($(it)): pot=$(rpot), kin=$(rkin)")
            else
                # <Nt| ρ |Perm(1)>
                _rpot, _rkin = spinful_rho_exch(_nupt, ndnt, rupt, rdnt, nebds1, volume;
                λe=λe, μee=μee, μep=μep, Cne=Cne)
                #println("ρ̃($(it)): pot=$(rpot), kin=$(rkin)")
                #rpota, rkina = spinful_rho_at(_nupt, ndnt, rupt, rdnt, nebds1, volume;
                #λe=λe, μee=μee, μep=μep, Cne=Cne)
                #println("ρ($(it)): pot=$(rpota), kin=$(rkina)")
            end
            #判断是否可以更新
            #if true
            # TODO: 这里要判断
            #ρ_old = exp(-E_old) -> ρ_new = exp(-E_new)
            # R = ρ_new / ρ_old = exp(- (E_new-E_old))
            ρ_rel = (_lkin + _lpot) - lkin - lpot
            ρ_rel += (_rkin + _rpot) - rkin - rpot
            if rand() < exp(-ρ_rel)
                #println("<$(it-1)|ρ|$(it)> = $(lpot), $(lkin) -> <$(it-1)|ρ|$(it)> = $(_lpot), $(_lkin)")
                #println("<$(it)|ρ|$(it+1)> = $(lpot), $(lkin) -> <$(it-1)|ρ|$(it)> = $(_lpot), $(_lkin)")
                nupt[ip, :] .= _nupt[ip, :]
                lpot = _lpot
                lkin = _lkin
                rpot = _rpot
                rkin = _rkin
            end
        end
        #把位置更新到WrdFb当中
        for ip in Base.OneTo(Ne)
            fbup[ip][it, :] .= nupt[ip, :]
        end
        nupt = beads(fbup, it)
        #println("beads after $(nupt)")#dn
        #尝试更新dn的位置
        #
        for ip in Base.OneTo(Ne)
            _ndnt = copy(ndnt)
            ran = rand(3).-0.5
            ran = ran*0.1
            #println(ran)
            _ndnt[ip, :] = map_to_1st(_ndnt[ip, :] + (ran)*λe, volume)
            #更新后的密度矩阵
            if it != 1
                #<t-1| ρ |t>
                _lpot, _lkin = spinful_rho_at(lupt, ldnt, nupt, _ndnt, nebds1, volume;
                λe=λe, μee=μee, μep=μep, Cne=Cne)
            else
                #<Nt| ρ |Perm(1)> 
                _lpot, _lkin = spinful_rho_exch(lupt, ldnt, nupt, _ndnt, nebds1, volume;
                λe=λe, μee=μee, μep=μep, Cne=Cne)
            end
            #右侧的密度矩阵
            if it != Nt
                # <t| ρ |t+1>
                _rpot, _rkin = spinful_rho_at(nupt, _ndnt, rupt, rdnt, nebds1, volume;
                λe=λe, μee=μee, μep=μep, Cne=Cne)
            else
                # <Nt| ρ |Perm(1)>
                _rpot, _rkin = spinful_rho_exch(nupt, _ndnt, rupt, rdnt, nebds1, volume;
                λe=λe, μee=μee, μep=μep, Cne=Cne)
            end
            #
            ρ_rel = (_lkin + _lpot) - lkin - lpot
            ρ_rel += (_rkin + _rpot) - rkin - rpot
            if rand() < exp(-ρ_rel)
                #println("<$(it-1)|ρ|$(it)> = $(lpot), $(lkin) -> <$(it-1)|ρ|$(it)> = $(_lpot), $(_lkin)")
                #println("<$(it)|ρ|$(it+1)> = $(lpot), $(lkin) -> <$(it-1)|ρ|$(it)> = $(_lpot), $(_lkin)")
                ndnt[ip, :] .= _ndnt[ip, :]
                lpot = _lpot
                lkin = _lkin
                rpot = _rpot
                rkin = _rkin
            end
        end
        #把位置更新到WrdFb当中
        for ip in Base.OneTo(Ne)
            fbdn[ip][it, :] .= ndnt[ip, :]
        end
        ndnt = beads(fbdn, it)
    end
    #更新结束
end


"""
计算当前的相位
"""
function abs_phase(fbup::WrdFb{T}, fbdn::WrdFb{T}, volume::Vector{Float64};
    λe=_ThermWL) where T
    Ne = length(fbup)#dn
    Nt, Nd = size(fbup[1])#dn
    #
    up0 = beads(fbup, 1)
    upNt = beads(fbup, Nt)
    emup = exch_free_kinetic(upNt, up0, volume; λe=λe)
    #
    dn0 = beads(fbdn, 1)
    dnNt = beads(fbdn, Nt)
    emdn = exch_free_kinetic(dnNt, dn0, volume; λe=λe)
    #
    detup = det(emup)
    detdn = det(emdn)
    sgn = detup*detdn
    pha = real(sgn / abs(sgn))
    return pha
end



"""
计算内能
"""
function internal_energy(fbup::WrdFb{T}, fbdn::WrdFb{T}, fbne::WrdFb{T},
    volume::Vector{Float64}, Δβ::Float64;
    λe=_ThermWL, μee=_WLElec, μep=_WLNeur, Cne=_NeurCoul) where {T}
    Ne = length(fbup)#dn
    Nt, Nd = size(fbup[1])#dn
    Nc = length(fbne)
    #对粒子数
    eng_pt = 1.5 * (Nc + Ne*2)
    #质子之间的势能
    eng_cc = 0.0
    cpos = beads(fbne, 1)
    for i1 in Base.OneTo(Nc); for i2 in Base.OneTo(Nc)
        if i1 == i2
            continue
        end
        qij = (cpos[i1, :] - cpos[i2, :])
        qij = map_to_1st(qij, volume)
        squqij = sum(qij.^2)
        eng_cc += 0.5*(Cne^2)/sqrt(squqij)
    end; end
    #电子之间
    eng_ee = 0.0
    uppos = beads(fbup, 1)
    dnpos = beads(fbdn, 1)
    e0pos = vcat(uppos, dnpos)
    for it in Base.OneTo(Nt)
        uppos = beads(fbup, it)
        dnpos = beads(fbdn, it)
        epos = vcat(uppos, dnpos)
        for i1 in Base.OneTo(2*Ne); for i2 in Base.OneTo(2*Ne)
            if i1 == i2
                continue
            end
            rij = epos[i1, :] - epos[i2, :]
            r0ij = e0pos[i1, :] - e0pos[i2, :]
            #squrij = sum(rij.^2)
            #eng_ee += 0.5/sqrt(squrij)
            eng_ee += 0.5*Kelbg_repulse_dbeta(rij, r0ij; λ11=μee)
        end; end
    end
    eng_ee = eng_ee / Nt
    #交换能
    upr1 = beads(fbup, Nt)
    upr2 = beads(fbup, 1)
    upex = exch_free_kinetic(upr1, upr2, volume; λe=λe)
    updm = exch_derv_beta(upr1, upr2, volume; λe=λe)
    invup = inv(upex)
    dlnupdb = tr(updm*invup)/Δβ
    #
    dnr1 = beads(fbdn, Nt)
    dnr2 = beads(fbdn, 1)
    dnex = exch_free_kinetic(dnr1, dnr2, volume; λe=λe)
    dndm = exch_derv_beta(dnr1, dnr2, volume; λe=λe)
    invdn = inv(dnex)
    dlndndb = tr(dndm*invdn)/Δβ
    eng_exch = dlnupdb + dlndndb
    #
    #电子和质子之间的势能
    cpos = beads(fbne, 1)
    uppos = beads(fbup, 1)
    dnpos = beads(fbdn, 1)
    e0pos = vcat(uppos, dnpos)
    eng_ce = 0.0
    for it in Base.OneTo(Nt)
        uppos = beads(fbup, it)
        dnpos = beads(fbdn, it)
        epos = vcat(uppos, dnpos)
        for ie in Base.OneTo(2*Ne); for ic in Base.OneTo(Nc)
            x0ij = cpos[ic, :] - e0pos[ie, :]
            xij = cpos[ic, :] - epos[ie, :]
            eng_ce += Kelbg_attract_dbeta(xij, x0ij; λ22=μep, Cne=Cne)
        end; end
    end
    eng_ce = eng_ce / Nt
    #βE = eng_pt + β*(eng_cc + eng_ee + eng_ce - eng_exch)
    return eng_cc, eng_ee, eng_ce, eng_exch
end

